/*
 * DoS.h
 *
 *  Created on: 2011-1-8
 *      Author: zhang
 */

#include <math.h>
#include <fstream>
#include <string>
#include <compiler_conf.h>

namespace KPM{
	//create dos graph data
	//template<class T>
	//void createDoS(T* h_U,int M,long hits,char* fileName);
	template<class T>
	class KPMDOS{
	private:
		double* e;
		T* d;
		unsigned int len;
		KPMDOS(){}
		/*******************************
		 * create kernel function, jackson kernel for default.
		 * M: number of moments.
		 ****************************/
		HOST STATIC INLINE T* createKernel(int M);

		/******************************
		 * this function create a the dos just for a single point.
		 * mu: chebyshev moments array
		 * m: number of chebyshev moments
		 * x: the point at which the KPMDOS value would be calculated.
		 ****************************/
		HOST STATIC T createAtSinglePoint(T* mu, unsigned int M, T* kernel, T x);

		HOST void restore(Scale<T>& scale);
	public:
		void writeToFile(string s);
		static T* createFromMoments(T* mu, int M, int s);
		static void free(KPMDOS<T> *o);
	};

	template <class T>
	void KPMDOS<T>::free(KPMDOS<T> *o){
		free(o->d);
		free(o->e);
		delete o;
	}

	template <class T>
	T* KPMDOS<T>::createKernel(int m){
		T* mu = (T*)malloc(sizeof(T) * m);
		for(int j=0;j<m;j++){
			mu[j] = (T)1/(m + 1)*((m-j+1)*cos(KPI*j/(m+1))
										   + sin(KPI*j/(m+1))*cos(KPI/(m+1))/sin(KPI/(m+1)));
		}
		return mu;
	}

	template <class T>
	T KPMDOS<T>::createAtSinglePoint(T* mu, unsigned int m, T* kernel, T x)
	{
		unsigned int j=0;
		T Re=0, T0, T1, T2;

		//j = 0;
		T0 = 1;
		Re += mu[0] * kernel[0];

		//j = 1;
		T1 = x;
		Re += 2 * mu[1] * T1 * kernel[1];

		for(j=2;j<m;j++){
		  T2 = 2 * x * T1 - T0;
		  Re += 2 * mu[j] * T2 * kernel[j];
		  T0=T1;
		  T1=T2;
		}
		Re *= 1/(3.14*sqrt(1-x*x));
		return Re;
	}


	/*****************************
	 *
	 * this function create dos function
	 * mu: chebyshev moments
	 * m: number of chebyshev moments
	 * s: number of samples selected for calculating KPMDOS
	 *
	 *****************************/
	template<class T>
	KPMDOS<T>* KPMDOS<T>::create(T* mu, unsigned int m, int s)
	{
		KPMDOS<T>* dos = new KPMDOS();
		dos->e = (T*)malloc(sizeof(T) * s);
		dos->d = (T*)malloc(sizeof(T) * s);
		memset(dos->d,0,sizeof(T) * s);

		T* kernel = KPMDOS<T>::createKernel<T>(m);

		for(int ei=0;ei<s;ei+=1){
			T x=(T) (ei-s/2) / (s+1) * 2;
			dos->d[ei]=KPMDOS<T>::createAtSinglePoint<T>(mu, m, kernel, x);
			dos->e[ei] = x;
		}
		free(kernel);
		dos->len = s;
		return dos;
	}

	template<class T>
	void KPMDOS<T>::writeToFile(string s){
		ofstream out(s);
		if(!out){
			for(int i=0;i<len;i++){
				out<<setpxecision(5)<<this->e[i]<<"\t"<<this->d[i]<<endl;
			}
		}
	}
}
